Spinodal amplification of density fluctuations 
in fluid-dynamical simulations of relativistic nuclear collisions 
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Extending a previously developed two-phase equation of state, we simulate head-on relativistic 
lead-lead collisions with fluid dynamics, augmented with a finite-range term, and study the effects 
of the phase structure on the evolution of the baryon density. For collision energies that bring the 
bulk of the system into the mechanically unstable spinodal region of the phase diagram, the density 
irregularities are being amplified significantly. The resulting density clumping may be exploited as 
a signal of the phase transition, possibly through an enhanced production of composite particles. 

PACS numbers: 25.75.-q, 47.75. +f, 64.75. Gh, 81.30.Dz 



Strongly interacting matter is expected to posses a rich 
phase structure. In particular, compressed baryonic mat- 
ter may exhibit a first-order phase transition that persists 
up to a certain critical temperature [l[ and experimental 
efforts are underway to search for evidence of this phase 
transition and the associated critical end point 0t3]- 

For this endeavor to be successful, it is important to 
identify observable effects that may serve as signals of the 
phase structure. This is a challenging task because the 
colliding system is relatively small, non-uniform, far from 
global equilibrium, and rapidly evolving, features that 
obscure the connection between experimental observables 
and the idealized uniform equilibrium matter described 
by the equation of state. Therefore, to understand how 
the presence of a phase transition may manifest itself 
in the experimental observables, it is necessary to carry 
out dynamical simulations of the collisions with suitable 
transport models. 

A first-order phase transition introduces spinodal in- 
stabilities [5-8] and the associated non-equilibrium dy- 
namics may generate observable fluctuations in the chiral 
order parameter 0, [l(| and the baryon density pdl - H^ ] . 
In order to explore this latter prospect, we simulate nu- 
clear collisions with finite-density fluid dynamics, using a 
previously developed two-phase equation of state. Fluid 
dynamics has the important advantage that the equation 
of state appears explicitly, contrary to most microscopic 
transport models where it is often unknown or is cumber- 
some to determine. Because it is essential to incorporate 
finite-range effects when seeking to describe the spinodal 
phase decomposition [1, H, H uM > we introduce a gradi- 
ent term into the expression for the local pressure in a 
non-uniform medium. This term ensures that two co- 
existing bulk phases will develop a diffuse interface and 
we calculate the associated temperature-dependent ten- 
sion. The gradient term also causes the dispersion rela- 
tion for the collective modes in the unstable phase region 
to exhibit a maximum, as is a characteristic feature of 
spinodal decomposition Q. Thus we obtain a transport 
model that has an explicitly known two-phase equation 
of state and that treats the associated physical instabil- 



ities in a numerically reliable manner. This is the first 
time a transport model with these key characteristics has 
been developed for high-energy nuclear collisions. 

As a first application of this novel tool, we simulate 
head-on collisions of lead nuclei at various energies. For 
a certain window of collision energies, several GeV per 
nucleon, the bulk of the system will reside within the 
spinodal region of the phase diagram for a sufficient time 
to allow the associated instabilities to enhance the initial 
density irregularities, a phenomenon that was exploited 
previously to obtain experimental evidence of the nuclear 
liquid-gas phase transition [Sllla]. As a quantitative mea- 
sure of the effect, we extract the moments of the density 
distribution and compare the degree of enhancement with 
what would result in the absence of the phase transition. 

In order to obtain a suitable equation of state, we em- 
ploy the method developed in Ref. [8j. Thus we work 
(at first) in the canonical framework and, for a given T, 
we obtain the free energy density /t(p) in the phase co- 
existence region by performing a suitable spline between 
two idealized systems (either a gas of pions and inter- 
acting nucleons or a bag of gluons and quarks) held at 
that temperature. In Ref. [8| the focus was restricted 
to subcritical temperatures, T < Tcrit, so for each T the 
spline points were adjusted so the resulting fr{p) would 
exhibit a concave anomaly, i.e. there would be two den- 
sities, pi(T) and p2(T), for which the tangent of fr(p) 
would be common. This ensures phase coexistence be- 
cause then the chemical potentials px = dpfrip) match, 
Pr(pi) = Mr (pa), as do the pressures, pr(pi) = Pripi)- 
In the present work, we have extended the equation of 
state to T > T cr it by using splines that have no concavity, 
as is characteristic of single-phase systems. After having 
thus constructed /t(p) for a sufficient range of T and p, 
we may obtain the pressure, pr{p) = pdpfr(p) — /t(p), 
and the energy density, s T (p) = fc(p) - Tdrfrip), b Y 
suitable interpolation and then tabulate the equation of 
state, po(s, p), on a convenient cartesian lattice. 

The resulting phase diagram in the (e, p) phase plane is 
illustrated in Fig.[TJ At a given net baryon density p, the 
lower bound on the energy density is given by £t=o(p). 
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Net baryon density p (in units of p ) 



FIG. 1: [Color online]The (e,p) phase diagram: St=o(p), the 
minimal e for a given p (black); the phase-coexistence bound- 
ary (red) ending at the critical point (c.p.), as well as the 
isothermal spinodal boundary where vt — (green) and the 
isentropic spinodal boundary where v 3 = (blue). Points 
corresponding to counter-streaming Lorentz-contracted nuclei 
are shown for various energies i?iab, indicated in ylGeV (cir- 
cles). The two contours (dashed) show where the points of 
maximum compression are concentrated in ensembles of col- 
lisions simulated with fluid dynamics for _Bi a b=2 or 5 AGeV. 



which represents the effect of compression. The lower 
and upper boundaries of the phase coexistence region 
are traced out by pi (T) and p 2 (T) , respectively. Uniform 
matter is thermodynamically unstable inside this region; 
while it is mechanically metastable near the phase co- 
existence boundaries, it looses mechanical stability when 
the speed of sound vanishes and density irregularities are 
amplified inside this spinodal boundary. 

It is important to recognize that the model describes 
the instabilities not only in the unstable spinodal region, 
but also those in the surrounding metastable region in 
which finite seeds are required for amplification to oc- 
cur (yielding nucleation or bubble formation jl3|). Thus 
density irregularities may be amplified by the metastable 
region as well, and any clumping generated inside the 
spinodal region may be further enhanced as the system 
expands through the nucleation region. 

When the dynamical evolution is governed by ideal 
fluid dynamics (see below) the instability boundary is 
characterized by the vanishing of the isentropic sound 
speed, v s = 0, where v 2 = (p/h) (dp/dp) s , p , with s be- 
ing the entropy density and h = p + e the enthalpy den- 
sity. For dissipative evolutions the finite heat conductiv- 
ity causes the instability region to widen as the boundary 
is then characterized by the vanishing of the isothermal 
sound speed, vt = 0, where v\ = (p/h) (dp/dp) T < v 2 . 
These instability boundaries are also indicated on Fig. [T] 

In order to ascertain the dynamical effects of the first- 
order phase transition, we construct a one-phase part- 
ner equation of state by replacing, for any T < T cr it, 



the actual fr(p) by the associated Maxwell construction, 
f¥(P) = fr(pi) + (p~ Pi)p(Pi), through the coexistence 
region Pl (T) <p< p 2 (T) where f^(p) < f T (p). 

For our present investigation, we describe the evolution 
of the colliding system by ideal fluid dynamics, because 
dissipative effects are not expected to play a decisive role 
for the spinodal clumping [15j . The basic equation of 
motion expresses four-momentum conservation, dfjT^ = 
0, where T^(x) = \p(x) + e(x)]u^(x)u u (x) - p(x)g» v 
is the stress tensor. It is supplemented by the continuity 
equation d^N^ = for the baryon charge current density 
N ,J, (x) = p(x)u fl (x), where u^(x) is the four- velocity of 
the fluid. These equations of motion are solved by means 
of the code SHASTA [17[ in which the propagation in the 
three spatial dimensions is carried out consecutively. 

A proper desciption of spinodal decomposition requires 
that finite-range effects be incorporated 0, Q . Therefore, 
following Rcfs. [§. Il5|. we write the local pressure as 

P(r) = Po(e(r), p(r)) - a 2 -| p(r)V 2 p(r) , (1) 

pi 

where we recall that p$(e, p) is the equation of state, the 
pressure in uniform matter characterized by e and p. Fur- 
thermore, p s = 0.153/fm 3 is the nuclear saturation den- 
sity and e s ~ m^p s is the associated energy density. 
The strength of the gradient term is then conveniently 
governed by the length parameter a. 

The gradient term will cause a diffuse interface to de- 
velop when matter of two coexisting phases are brought 
into physical contact. The associated interface tension 
can readily be determined from the equation of state Q , 

r 2(T) mi 

u T = a / [2e s Af T (c)] 1/2 dc, (2) 

J Cl (T) 

where c = p/ p s denotes the degree of compression and 
Af T (c) = /r(c) — fri ) IS the difference between the 
actual free energy density at the specified compression 
and that associated with the corresponding Maxwell con- 
struction. The tension ot decreases steadily as T is 
raised and finally vanishes at T cr it; it is displayed in Fig. 
[5] as obtained for various values of the range a. 

Furthermore, uniform matter inside the spinodal re- 
gion (where v 2 < 0) is mechanically unstable and den- 
sity ripples of wave number k will be amplified at a rate 
7fc(p, e) that depends on the strength of the gradient 
term, y% — |u s | 2 fc 2 — a 2 (e s /h)(p/ p s ) 2 k 4 . Thus the gra- 
dient term introduces a penalty for the development of 
short-range undulations and 7fc exhibits a maximum at 
the favored length scale; the familiar dispersion relation 
for ideal fluid dynamics, 7^ = \v s \k, is recovered in the 
absence of a gradient term, a = 0. 

These spinodal growth rates have been extracted by 
following the fluid-dynamical evolution of small harmonic 
perturbations in uniform matter. Figure [2] shows the re- 
sulting 7^ curves extracted at the phase point (6p s , 10e s ) 
located in the central spinodal region (see Fig. Q]). Be- 
cause the numerical algorithm inevitably produces some 
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FIG. 2: [Color online] Effect of the gradient term: For various 
values of the parameter a in Eq. (JTJ is shown the interface 
tension ctt as a function of T (upper panel) and the growth 
rate 7fc as a function of the wave number k of a disturbance, 
as obtained with ideal fluid dynamics (lower panel). 



degree of dissipation [17| , the growth rates extracted for 
a = deviate from the exact ideal dispersion relation 
7fe = \v s \k. Our studies suggest that while those results 
are reliable only up to k < 4/fm the use of a finite range, 
a > 0, extends validty to significantly finer scales. 

Having thus shown that our model produces both 
a meaningful interface tension and reasonable spinodal 
growth rates, we apply it to central collisions of lead 
nuclei bombarded onto a stationary lead target at var- 
ious kinetic energies, -Ef a b- For each energy, an ensemble 
of 200 separate evolutions are generated, each starting 
from a different initial state obtained by means of the 
UrQMD model (in its cascade mode) [18| - |20| in order to 
take account of the non-equilibrium dynamics in the very 
early stage of the collision. The switch from UrQMD to 
fluid dynamics occurs at a time to when the two Lorentz- 
contracted nuclei have just passed through each other 
and all initial collisions have occurred. The needed fluid- 
dynamical quantities, p(r), e(r), u(r), are then extracted 
from the UrQMD state by representing each hadron by 
a gaussian of 1 fm width and mapping this information 
onto a cartesian spatial lattice [21] . Thus the effects of 



FIG. 3: [Color online] Mean time evolution of normalized den- 
sity moments obtained for E^b = 3A GeV with a = 0.033 fm 
(for which ao ~ 10 MeV/fm 2 ), using either the two-phase 
equation of state (solid) or its one-phase partner (dashed). 

stopping as well as local event- by- event fluctuations are 
naturally included in the ensemble of the initial states. 

Figure [2] brings out the intimate relationship between 
the interface tension <tt and the spinodal growth rates 7^ 
and the range parameter a merely presents a convenient 
means for linking the two properties. For the present 
paper, we have adjusted a so that the zero-temperature 
interface tension is approximately 10 MeV/fm 2 , the value 
that was also preferred in Ref. [lj]. This corresponds to 
the middle curves in Fig. [5J obtained for a — 0.033 fm. 

The key effect of the first-order phase transition is the 
spinodal amplification of spatial irregularities [!, 0, H, 
Il5| . In order to obtain a quantitative measure of the 
resulting degree of clumping in the system, we extract the 
moments of the (net) baryon density distribution p(r), 

(p N ) ee -L J p ( r ) N p(r)d 3 r, (3) 

where A = j p(r)d 3 r is the total (net) baryon number. 
Thus the corresponding normalized moment, (p N )/ (p) N , 
is unity for N — 1. These quantities have observational 
relevance because they are intimately related to the rel- 
ative production yield of composite baryons. 

The time evolution of the normalized density moments 
is illustrated in Fig. [31 Generally, for a given density 
distribution p(r), the normalized moments increase with 
the order TV, as one would expect. At the initial time 
to the pre-equilibrium fluctuations are transcribed into 
the fluid-dynamical functions (see above) and they may 
then be amplified subsequently by the instabilities as- 
sociated with the first-order phase transition. For com- 
parison we show also the corresponding results for the 
one-phase partner equation of state that has no insta- 
bilities. There is a striking difference between the two 
sets of results: Whereas the one-phase equation of state 
hardly produces any amplification at all, the one having 
a phase transition leads to significant enhancements that 
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FIG. 4: [Color online] Mean maximum enhancement of the 
normalized density moments for N = 7,5 (top, bottom) as ob- 
tained for various energies using either the two-phase equation 
of state (solid) or its one-phase Maxwell partner (dashed). 

increase progressively with N and amount to well over an 
order of magnitude for N = 7. This clearly demonstrates 
that the first-order phase transition may have a qualita- 
tive influence on the dynamical evolution of the density. 
Due to the variations in the initial conditions, the phase 
regions explored differ from one collision to the other. As 
a result, some of the evolutions experience amplifications 
that are significantly larger than the ensemble average 
(by up to a factor of five or so) , whereas more evolutions 
are affected considerably less. 

Figure [3] depicts what would happen if the expansion 
were to continue while maintaining local equilibrium so 
the generated density enhancements eventually subside. 
However, as the hadronic gas grows more dilute, local 
equilibrium cannot be maintained. Consequently, if the 
decoupling occurs sufficiently soon after the clumps are 
formed, the associated phase-space correlations may sur- 
vive. Studies of this issue are underway. 

The degree of density clumping generated during a col- 
lision depends on how long time the bulk of the matter 
is exposed to the spinodal instabilities. The optimal sit- 
uation occurs for beam energies that produce maximum 
bulk compressions lying well inside the unstable phase re- 
gion because the instabilities may then act for the longest 
time [1, UU- At lower energies an ever smaller part of 
the system reaches instability and the resulting enhance- 
ments are smaller. Conversely, at higher energiers the 
maximum compression occurs beyond the spinodal phase 
region and the system is exposed to the instabilities only 
during a relatively brief period during the subsequent ex- 
pansion. For still higher energies the spinodal region is 
being missed entirely. 

Figure 2] shows the (ensemble average) maximum en- 
hancement achieved as a function of the beam energy 



for the two equations of state. The existence of an op- 
timal collision energy is clearly brought out. While the 
presently employed equation of state suggests that this 
optimal range is E\ a ^ ss 2 — 4 A GeV, it should be rec- 
ognized that others may lead to different results. Fur- 
thermore, the magnitude of the effect depends on the de- 
gree of fluctuation in the initial conditions which in turn 
are governed by the pre-equilibrium dynamics. On the 
other hand, our studies suggest that the optimal energy 
is rather insensitive to the range parameter a. 

To summarize, by augmenting an existing fluid- 
dynamical code with a gradient term, we have obtained 
a transport model that is suitable for simulating nuclear 
collisions in the presence of a first-order phase transition: 
it describes both the tension between coexisting phases 
and the dynamics of the unstable spinodal modes. Apply- 
ing this novel model to lead-lead collisions, we have found 
that the associated instabilities may cause significant am- 
plification of initial density irregularities. Because such 
clumping is expected to facilitate the formation of com- 
posite particles, such as deuterons and tritons, this effect 
may be experimentally observable. 

However, the magnitude of the generated clumping de- 
pends considerably on the largely unknown specifics of 
the equation of state which determines whether the un- 
stable region of the phase diagram is entered during the 
collision, for how long the system remains there, and how 
rapidly irregularities are amplified. It would therefore be 
interesting to expand the present kind of study to other 
equations of state and also to ascertain the importance 
of dissipative effects [H, HH . 

Furthermore, standard fluid dynamics propagates the 
system deterministically and thus ignores the effect of 
the ever present thermal noise. Because the initial state, 
in the present study, is already endowed with significant 
fluctuations, the thermal noise is likely to play only a 
minor role (see Ref. Q), but it would be interesting to 
investigate this quantitatively. (A study of the effect of 
noise in chiral fluid dynamics was recently made [Hj].) 

Spinodal decomposition is an inherent characteristic 
of first-order phase transitions and this non-equilibrium 
phenomenon may therefore provide especially com- 
pelling signals. The present study demonstrates that the 
phase structure does affect the character of the density 
evolution and we hope that these results will stimulate 
efforts to develop analysis techniques for extracting the 
related observables from experimental data. 

Acknowledgments 

We wish to acknowledge stimulating discussion with 
Volker Koch. This work was supported by the Office of 
Nuclear Physics in the U.S. Department of Energy's Of- 
fice of Science under Contract No. DE-AC02-05CH11231; 
JS was supported in part by the Alexander von Humboldt 
Foundation as a Feodor Lynen Fellow. 



5 



[1] M.A. Stephanov, K. Rajagopal, and E.V. Shuryak, Phys. 

Rev. Lett. 81, 4816 (1998). 
[2] B.I. Abelev et al, Internal STAR Note SN0493, 2009: 

drupal.star.bnl.gov /STAR/starnotes / public /sn0493. 
[3] B. Friman, C. Holme, J. Knoll, S. Leupold, J. Ran- 

drup, R. Rapp, and P. Senger, The CBM Physics Book, 

(Springer, Berlin, 2010). 
[4] A.N. Sissakian and A.S. Sorin, J. Phys. G 36, 064069 

(2009). 

[5] Ph. Chomaz, M. Colonna, and J. Randrup, Phys. Rep. 

389, 263 (2004). 
[6] J. Randrup, Phys. Rev. Lett. 92, 122301 (2004). 
[7] C. Sasaki, B. Friman, and K. Redlich, Phys. Rev. Lett. 

99, 232301 (2007). 
[8] J. Randrup, Phys. Rev. C 79, 054911 (2009). 
[9] K. Paech, H. Stocker, and A. Dumitru, Phys. Rev. C 68, 

044907 (2003). 

[10] M. Nahrgang, S. Leupold, C. Herold, and M. Bleicher, 
Phys. Rev. C 84, 024912 (2011). 



[11] J. Randrup, Acta Phys. Hung. A 22, 69 (2005). 

[12] V. Koch, A. Majumder, and J. Randrup, Phys. Rev. C 

72, 064903 (2005). 
[13] I.N. Mishustin, Phys. Rev. Lett. 82, 4779 (1999). 
[14] D. Bower and S. Gavin, Phys. Rev. C 64, 051902 (2001) 
[15] J. Randrup, Phys. Rev. C 82, 034902 (2010). 
[16] B. Borderie et al, Phys. Rev. Lett. 86, 3252 (2001). 
[17] D.H. Rischke, S. Bernard, and J. A. Maruhn, Nucl. Phys. 

A 595, 346 (1995). 
[18] S.A. Bass et al, Prog. Part. Nucl. Phys. 41, 225 (1998). 
[19] M. Bleicher et al, J. Phys. G 25, 1859 (1999). 
[20] H. Petersen, J. Steinheimer, G. Burau, M. Bleicher, and 

H. Stocker, Phys. Rev. C 78, 044901 (2008). 
[21] J. Steinheimer et al, Phys. Rev. C 77, 034901 (2008). 
[22] V.V. Skokov and D.N. Voskresensky, Nucl. Phys. A 828, 

401 (2009). 

[23] M. Nahrgang, C. Herold, S. Leupold, I. Mishustin, and 
M. Bleicher. larXiv:1105.1962l [nucl-th]. 



